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Abstract 

A computational study of morphological instabilities of a two-dimensional nematic front under 
directional growth was performed using a Landau-de Gennes type quadrupolar tensor order pa- 
rameter model for the first-order isotropic/nematic transition of 5CB (pentyl-cyanobiphenyl). A 
previously derived energy balance, taking anisotropy into account, was utilized to account for la- 
tent heat and an imposed morphological gradient in the time-dependent model. Simulations were 
performed using an initially homeotropic isotropic/nematic interface. Thermal instabilities in both 
the linear and non-linear regimes were observed and compared to past experimental and theoretical 
observations. A sharp-interface model for the study of linear morphological instabilities, taking 
into account additional complexity resulting from liquid crystalline order, was derived. Results 
from the sharp-interface model were compared to those from full two-dimensional simulation iden- 
tifying the specific limitations of simplified sharp-interface models for this liquid crystal system. In 
the non-linear regime, secondary instabilities were observed to result in the formation of defects, 
interfacial heterogeneities, and bulk texture dynamics. 
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I. INTRODUCTION 



Directional growth of materials undergoing a phase transition consists of pulling the 
material through a temperature gradient, usually from a higher to a lower temperature 
(relative to the phase transition temperature), such that an interface between the two phases 
is established in the central region of the sample. A schematic of this setup is shown in Fig. [IJ 
Directional growth configurations have practical uses such as producing single crystals with 
defined crystal orientation, producing eutectic composite materials, and as a purification 
method. In addition to these applications, directional growth configurations are also widely 
used to study the fundamental nature of phase transitions. 




Figure 1: Schematic of a directional growth experimental system where a slide/thin film enclosing 
a material is pulled at a constant velocity v from a plate or oven at a temperature above the 
phase transition temperature to one that is below the phase transition temperature; the opaque 
portion of the slide represents the ordered/solid phase and the translucent portion represents the 
disordered /liquid. 



One important characteristic present in general growth processes that involve phase tran- 
sitions, and in the directional growth system in particular, is the occurrence of morphological 
instabilities [1]. These instabilities arise in a great variety of materials undergoing different 
types of transitions, ranging from diffusive phase transitions to phase ordering. The study 



of these instabi 

n 



ities has seen substantial advances since the seminal work of Mullins and 
Sekerka 0, [j, |4j. Since then, some examples of systems that have been studied, in addi- 
tion to the solidification of binary afioys flflBQB, incfiade sificon wafer production 
in the semiconductor industry [9], polymer crystallization [lfj, and liquid crystal growth 



HQ, 
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The fundamental understanding of morphological instabilities is important to the science 
of phase transformations. These instabilities also affect the final material morphology (cel- 
lular, dendritic, textures in liquid crystals et al) and, subsequently, the structural and func- 
tional material properties. For example, the internal structure of a liquid crystal spherulite 
controls the optical performance of polymer-dispersed liquid crystals. Previous work jlol . 1(| 
has shown that liquid crystal spherulite growth on the nanoscale leads to a series of orien- 
tational events, including defect nucleation and shape instabilities, showing that a better 
understanding of growth laws may be useful to develop improved optical materials. 

The study of directional growth in liquid crystalline material systems has been the fo- 
cus of much previous work. These materials offer accessible time-scales for experimental 
observations as well as the combination of soft and anisotropic behavior. The anisotropic 
character of liquid crystals results in many secondary instabilities, in addition to those de- 
scribed by Mullins and Sekerka 

as, a. 

While liquid crystalline systems are templates 
for the overarching study of morphological instabilities, they are also pervasively used in 
industrial applications. Their typical use in thin-film geometries, where large temperature 
gradients can be easily produced, places further importance on the need to characterize and 
understand the occurrence of morphological instabilities. 

nnnn 



Much past experimental 
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14], [171, [18|, 119J and theoretical [12|, |20| work has focused 



on studies of thermal and other types of instabilities of isotropic/nematic mesophase tran- 
sitions in directional growth experiments. This system is a convenient starting point since 
the nematic phase is the simplest of the liquid crystal mesophases. Instabilities in smectic 
and columnar liquid crystals have been less studied 2l|]. Much of the work on morphologi- 
cal instabilities in nematic fronts has focused on explaining discre pan cies between capillary 



lengths, determined experimentally 
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191 ] and theoretically 



14 



tigations of the effects of convection and impurity concentration 



17 



22j. Experimental inves- 
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have been conducted, 



23| resulting from the 



but the diverse set of possible morphological instabilities 
inherent anisotropy and anchoring effects have not yet been completely explored. The ex- 
perimental study of these instabilities is inherently difficult due to the length and time scales 
involved. On the other hand, computational studies are able to access these scales and shed 
light on the governing principles. 

The overall mechanism which typically results in morphological instabilities is rooted in 
the existence of a temperature gradient, where the velocity of the interface v is proportional 
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to [15|]: 

v oc \AF\-C 

\AF\ocT s -T (1) 

where \AF\ is the free energy difference between the ordered/disordered phases, C is the 
capillary force, T s is the coexistence temperature of the two phases, and T is the temperature 
of the material. 

When a shape perturbation arises in a moving interface, in the presence of a temperature 
gradient (refer to Fig. |2J), it finds a relatively increased \AF\. Subsequently, this promotes 
perturbation growth, increasing its relative velocity Eq. (CQ) (with respect to the unperturbed 
front). In addition to this effect, the perturbation has curvature and thus a capillary force is 
exerted upon it, inhibiting growth. This thermo-capillary competition results in the growth 
of perturbations with a wavelength (or radius of curvature) greater than a critical value, 
which decreases as the temperature gradient increases. 



Solid f Liquid 

Figure 2: Schematic representation of the Mullins-Sekerka instability [2. !] where the bold 
vertical line represents a flat ordered/solid front growing into the disordered/liquid phase under 
an imposed temperature gradient, where the temperature profile is shown by the stippled line. A 
perturbance in the front growth, the dotted line, is both inhibited by capillary forces and promoted 
by a greater decrease in free energy from phase transition as temperature decreases (see Eq. (HJ); 
note that the temperature of the ordered phase, in the vicinity of the interface (as shown), is 
greater than in the disordered phase due to latent heat effects. 

Although not studied in this work, it is also possible to observe instabilities induced by 
gradients in concentration. This is referred to as constitutional super-cooling J24], and is 
the case in mixtures or when impurities are present. The ordered/solid phase has a lower 
equilibrium concentration of the solute or impurities than the disordered/liquid phase, thus 



the impurities are rejected to the disordered/liquid phase. This results in the concentration 
of the impurity at the interface being higher than in the bulk disordered/liquid phase. Any 
perturbation in the interface will find a lower concentration of impurity, resulting in a higher 



coexistence temperature and increased front velocity (re: 



24 



er to 
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jlq. CD). For a more in-depth 
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explanation of morphological instabilities see Refs. 

In a directional growth system, an appropriate temperature gradient, in the vicinity of 
the interface (see Fig. [2]), can be produced via latent heat resulting from the phase transi- 
tion. Thus, even though the externally imposed temperature gradient stabilizes the interface 
(higher temperature in the disordered/liquid phase than in the ordered/solid phase), a desta- 
bilizing temperature gradient can exist locally in the vicinity of the interface (depending on 
the magnitude of the heat of phase transition). For some value of the externally imposed 
temperature gradient, the local gradient at the interface will be exactly zero, thus for ex- 
ternal gradients greater than this value, the shape instability will not be observed. If the 
external gradient is lower than this value (or if it is negative), there will be some wavelength 
range for which the perturbations grow, so a shape instability is observed. At relatively low 
pulling velocity, stationary sinusoidal shape patterns are observed, produced by this tem- 
perature gradient. As velocity is increased, the patterns remain periodic but lose sinusoidal 
shape. At further increased velocities, nonlinear effects, including non-periodic/chaotic in- 



stabilities can be observed in some cases 



22j (see Fig. [3^). Finally, for very high velocities, 



a restabilization is possible which reforms a flat interface. 

Mullins and Sekerka were the first to model morphological instability using a sharp- 
interface model, in the linear regime, considering that growth was limited by diffusion 2], y]. 
In these models, the equations for heat or mass diffusion are solved in each phase, and 
the boundary is discontinuous. Utilizing this formulation, dispersion diagrams, where the 
growth velocity versus wavelength are plotted, are able to be obtained analytically for a 
sinusoidal perturbation. This type of model has been extended to account for some types 
of nonlinear phenomena 20, |22j, [28] and other effects, but are generally not feasible for use 
to model complex scenarios where nonlinear instabilities occur. 



Phase field models 



29| have also been used to model morphological instabilities. These 



models inherently capture a more complete set of physics involved in that the material is 
modeled as a continuum. Thus the interface is no longer assumed discontinuous and, in the 
case of liquid crystals, texturing processes can be resolved. Furthermore, interfacial and bulk 
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properties are determined implicitly by the model and its parameters. A general comparison 
of both types of models (for a scalar phase field model) can be found in Ref. ]&■ 

The Landau-de Gennes tensorial model for the isotropic/nematic transition 30] is one of 
the most effective theoretical approaches to capture the kinetics and dynamics of this tran- 



sition at mesoscopic scales 



31 



321 ]. The phenomenological nature of this model is conducive 



to experimental validation and, as a result of these comparisons, this model has been shown 
to be relatively successful at capturing the physics of the isotropic/nematic transition 31]. 
The more coarse-grained approach of the Landau-type models, allows access to multi-scale 
phenomena which are relevant to experimental observations, but at a resolution unattain- 
able other than through numerical simulation. Thus, this model has been applied to study 
a broad range of phase-ordering phenomena, from use as a template to study the formation 
of the early universe 33j to the structure of liquid crystalline fibers, membranes, films, and 



drops 34]. 

With respect to this model of isotropic/nematic liquid crystalline transition, Fig. [3b 
elucidates the key physics captured in a schematic of morphological instabilities in directional 
growth of nematic phase ordering; Q is the quadrupolar tensor order parameter, k is the 
outward unit normal, and b is the curvature tensor of the interface. In nematic liquid 
crystals, the interface shape is coupled to the temperature and the order parameter Q. The 
coupling between the degree of phase ordering and temperature is implicitly accounted for 
in the Landau-de Gennes tensorial model. The incorporation of a previously derived energy 



balance 



351 ] accounts for the heat of transition and anisotropy in heat conduction (arising 



from the imposed boundary temperatures and the heat of phase ordering). The Q-tensor 
model implicitly incorporates the non-trivial couplings between growth, shape, and texturing 
dynamics. 

The majority of theoretical approaches to the study of morphological instabilities in 
directional growth, both in the isotropic/nematic mesophase transition and the general 
case, utilize simple sharp-interface 22] models and focus mainly on the linear regime. The 
main advantage to this approach is that simple analytical solutions can be obtained and 
there is good agreement between the predictions of these types of models and experimental 
observations in the linear regime. These models are not suitable for the study of instabilities 
in the non-linear regime, where much of the physics that is neglected in linearized sharp- 
interface models plays a role. An example of these complex physics are, in the case of liquid 
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Figure 3: a) schematics of the different regimes of morphological instabilities where the control 
variables include the pulling velocity v and the imposed global temperature gradient AT//, where 
I is the distance between the heat sources b) schematic of the liquid crystal physics taken into 



351 ] for the isotropic/nematic 



account in the non-isothermal tensorial Landau-de Gennes model 
transition employed in this work: non-isothermal conditions (T1/T2 are the imposed boundary 
temperatures), bulk texture (Q B is the tensor order parameter field in the bulk), interfacial gradi- 
ents and heterogeneities (Q s is the tensor order parameter field at the interface), heat of transition 
(L), and anisotropy/capillary forces (k/b is the interface normal/curvature tensor). 



crystals, de 
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22 



ect formation and texturing processes, which have been experimentally observed 



23) . Thus, past theoretical work was focused on modeling the linear regime and 



some simple non-linear effects, but the complex texturing and defect formation in directional 
growth of liquid crystals has not been modelled before. The state of the art in the study of 
this seminal problems is summarized in Chapter B.VI of ref. 3(| where it is explained in 
detail that current simulation and theoretical work is limited to low-dimensional, director- 
type macroscopic models. In the present case multi-scaling and multi-dimensionality are 
extended to the current computational limits. 

The main objectives of this work on modeling the directional solidification of calamitic 
low molar mass nematic liquid crystals are: 

• compare the Q(x, y,t) and T(x,y,t) predictions of two-dimensional (x,y) transient 



simulation of the non-isothermal Landau-de Gennes model 



15 



30 



35] results in the lin- 



ear (small amplitude) regime with those from the standard analytical (sharp-interface) 



method. 



• use the non-isothermal Landau-de Gennes model to simulate Q(x, y,t) and T(x,y,t) 
in the nonlinear regime, where secondary instabilities result in texturing and defect 
formation. 

The main assumptions of this work are: 

• the physical set-up and geometry correspond to the classical directional solidification 
experiment, as shown in Fig. [TJ 

• the liquid crystal is a pure thermotropic calamitic low molar mass material (see Table 

• directional growth is considered to be driven by an externally imposed temperature 
gradient. 

• thermal fluctuations are neglected. 

• constant physical properties corresponding to the 5CB (pentyl-cyanobiphenyl) liquid 
16], including thermal conductivity js7], are used. 



crystal 



a periodic perturbation to the interface is assumed in the initial conditions. 



A full tensorial Landau-de Gennes model is used 



15 



30t ] and with an applicable previ- 



ously derived thermal energy balance 35( . The nonlinear sharp-interface theory is not used 
because analytical solutions are infeasible when taking into account anisotropy, defect for- 
mation, and textures. 

This work is organized in three distinct parts introducing the relevant physics of nematic 
directional growth, modeling and simulation using the tensorial Landau-de Gennes model, 
and presentation of the results. Sec. Ill Al presents a description of the isotropic/nematic 
phase transition and quadrupolar tensor order parameter used in the Landau-de Gennes 
model. Sec. Ill Bl presents a brief description of the Landau-de Gennes model and the 
nematic thermal energy balance (the latter accounts for dissipation, anisotropic heat con- 
duction, and latent heat of transition). Sec. Ill CI describes the simulation conditions and 
complexities resulting in the use of the high-order Landau-de Gennes model. Sec. [TTT] is 
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divided into two subsections presenting results of both linear and nonlinear regime mor- 
phological instabilities. The linear regime results were found through analytical solutions 
of a derived sharp-interface model and through simulation of the higher-order Landau de 
Gennes model. Nonlinear regime results were only accessible through simulation. Finally, 
Sec. HV] summarizes the findings and presents conclusions. These conclusions focus on the 
limitations to sharp-interface model approaches compared to the Landau-de Gennes model. 
Additionally, interesting results in the nonlinear regime are presented where a an interfacial 
reorientation is mechanism is found to occur through defect shedding. 

II. BACKGROUND AND THEORY 

A. Liquid crystalline order 

Liquid crystalline phases or mesophases are materials which exhibit partial orientational 
and/or translational order. They are composed of anisotropic molecules which can be disc- 
like (discotic) or rod-like (calamitic) in shape. Thermotropic liquid crystals are compounds 
that exhibit mesophase ordering in response to temperature changes. Lyotropic liquid crys- 
tals that most greatly exhibit mesophase behavior in response to concentration changes. 
Effects of pressure and external fields also influence mesophase behavior. This work focuses 
the study of calamitic thermotropic liquid crystals which exhibit a first-order mesophase 
transition. 

An unordered liquid, where there is neither orientational nor translational order (apart 
from an average intermolecular separation distance) of the molecules, is referred to as 
isotropic. Uniaxial nematic liquid crystalline order in rod-like mesogens involves partial 
orientational order and positional disorder, where an average orientational axis, known as 
the director, is observed. Schematic representations of these uniaxial nematic and isotropic 
ordering are shown in Fig. HI 

B. Order parameters and the phenomenological model 

Theoretical characterization of nematic order is accomplished using an order parameter 
that adequately captures the physics involved. This order parameter has an amplitude and 
phase associated with it. In order to characterize the partial orientational order of the 
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(a) 



(b) 



Figure 4: (Color online) schematics of the a) isotropic and b) uniaxial nematic phases with rod-like 
mesogens. 



nematic phase, a second order symmetric traceless tensor can be used 

Q = S (nn - -i) + (mm - 11) 



30] 



(2) 



where n/m/1 are the eigenvectors of Q, which characterize the average molecular orienta- 
tional axes, and S/P are scalars which represent the extent to which the molecules conform 



to the average orientational axes 



34 



38, 



391 ] . Uniaxial order is characterized by 5* and 



n, which correspond to the maximum eigenvalue (and its corresponding eigenvector) of Q, 
S = \nn- Biaxial order is characterized by P and m/l, which correspond to the lesser 
eigenvalues and eigenvectors, P — — | (/i m — fii). 



The Landau-de Gennes model for the first order isotropic/nematic phase transition 301 ] : 
f-fo = \a(Q : Q) ~\b{Q ■ Q) : Q + \c(Q : Q) 2 



+h 1 (VQ-.VQ) + h 2 (V-Q)-(V-Q) 
+\hQ ■ (VQ : VQ) 



(3) 



a = a (T-T NI ) (4) 

(5) 

where a, b, c are bulk parameters, T/v/ is the lower stability limit of the isotropic phase, 
and li, Z 2 , ^3 are the elastics constants. An equi-bend/splay assumption is used in this work, 

10 



resulting in h, 1% > and Z3 = 0. This assumption is made based upon previous studies of 
nematic spherulite morphology resulting from the interplay of splay, twist, and bend elastic 



constants 



16 



40 



4l) . Interfacial contributions implicitly result from the inclusion of the 
. Detailed past work has studied the interfacial contributions of 



gradients terms in Eq. 

ihese gradient terms including anchoring, curvature, and heterogeneous effects. See to Refs. 



15 



16 



40 



41] for a comprehensive study of these effects. 



The Landau-Ginzburg time dependent formulation [42] is used to minimize the free energy 
functional Eq. ([6]) of the simulation volume: 



fdV 



/'■ 



dQ_ 

dt 



8F_ 
6Q 



-1 ST 



(6) 
(7) 



where F is the total free energy, fi is the rotational viscosity and only the symmetric- 
traceless component of the functional derivative is utilized (denoted by the superscript ST). 
The general differential energy balance, neglecting convection, is [351 ] : 

dT dQ dQ f d df\ ST dQ 

pCp m =fi ^ : ^ + T \dQdfS : ^- v - q (8) 

where C p is the specific heat and q is the total heat flux. The first right hand side term 
in Eq. [8] is dissipation due to nematic order dynamics, the second the heat of transition 
to nematic order, and the last thermal diffusion. Temperature fluctuations are neglected in 
this model but could be incorporated via stochastic terms. 

The heat flux can be calculated from the anisotropic Fourier's law: 



K ■ VT 



(9) 



where the thermal conductivity tensor K is used due to the anisotropy of the nematic phase. 
The thermal conductivity tensor can be written as the sum of isotropic and anisotropic 
contributions: 

'k\\ + 2k±' 



K 



k iso d + k an Q = ( a " + 3 2A:± ) d + (k\\ - k ± ) Q 



(10) 



where ki. so and k an are the isotropic and anisotropic contributions to the thermal conduc- 
tivities, and k\\ and k± are the conductivities in the directions parallel and perpendicular to 
the director, respectively. 
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C. Simulation method 



A schematic of the geometry of the two-dimensional simulation domain and boundary 
condition types are shown in Fig. A central sub-domain, with a very refined mesh, is used 
in order to resolve the details at the nematic/isotropic interface. Two outer sub-domains, 
with coarser meshes, are used to resolve the gradients in the order parameter in the bulk 
nematic phase and in temperature. These bulk gradients are at length scales orders of 
magnitude greater than those at the interface. Thus, the use of the fine mesh in these sub- 
domains is not necessary to resolve the textures. In the right sub-domain, a bulk isotropic 
phase is assumed (verified a posteriori), only the thermal energy balance Eq. (jSJ) was solved. 

Due to computational limitations the length scales accessible via simulation in the present 
work are on the order of microns. Experimental observations of morphological instabilities 



111 . Il3l . Il4j . 1 171 Il8l . Il9| have observed characteristics length scales on the order of tens and 
hundreds of microns. To circumvent this, a destabilizing temperature gradient is imposed 
with both temperature boundary conditions (see Fig. [U boundary conditions 1 and 4) 
below the bulk transition temperature for the material parameters used. This allows for the 
observation of morphological instabilities at wavelengths accessible in the domain sizes. 

For symmetry boundary conditions of Q, vector symmetry considerations result in the 
the following invariants [351 ]: 

OQxx dQyy Q .Q = Q = Q.^£ = Q 

where Xi is the coordinate associated with the basis vector normal to the symmetry axis. 

A reference system moving with velocity v, in the x direction and equal to the negative 
of the velocity of the moving isotropic/nematic interface, was used. This maintains a static 
interface position, where time derivatives are replaced by material derivatives: 

d d d . . 

Initial conditions resulting in an initially flat unperturbed homeotropic interface were 
achieved using the Heaviside step function: 

Q = Q He(y-y ) (13) 

where Qo is the bulk Q-tensor value corresponding to the interface temperature, He() is the 
Heaviside step function, and y is the initial position of the interface. Simulations were run 
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Figure 5: Schematic of the computational domain where the left sub-domain is the bulk homoge- 
neously oriented bulk nematic domain, the center sub-domain encompasses the isotropic/nematic 
interface, the right sub-domain is fully isotropic, and the width/height of the domain are 2B and 
Z, respectively; numbering corresponds to the boundary conditions used where: 1) T = T\ and 
Neumann conditions are used for Q, 2) T/Q are coupled 3) T is coupled and Q = 0, 4) T = T2 and 
Neumann conditions are used for Q, and 5) symmetry boundary conditions for T/Q, see Eq. (jlip : 
note that in this work both T\ and T2 are below the bulk isotropic/nematic transition temperature 
and T2 < T\ due to computational limitations (see Sec. Ill C|) . 

to determine steady-state, or stationary, interface profiles. Using these computed profiles a 
perturbation condition was applied using the following transformation: 

Q(x) = Qstat [y + a cos (j-x^J e~ 5y \ (14) 

where Q s tat(v) is the stationary (computed) profile of Q, a is the initial amplitude of the 
sinusoidal perturbation, A the initial half-wavelength, and the coefficient 6 is the decay 
factor. The characteristic decay length, <5 -1 , was chosen to equal 200nm so that the pertur- 
bation decays rapidly in the bulk relative to initial length scale of the nematic domain (on 
the order of microns). 

The initial condition Eq. ( |T4l) represents a flat stationary interface with a small sinusoidal 
perturbation. A half wavelength was resolved in the simulation domain, where symmetry 
was used to capture a periodicity and the simulation domain size I was varied, to represent 
different wavelengths. Due to the limitations of the finite-element software used (Comsol 

13 



BC 5 





Isotropic 




Domain 




Interface 




Bulk Nematic 




Domain 



Multiphysics), a static mesh was iteratively determined for each simulation where mesh 
density ranged from a maximum of approximately one second-order Lagrange element per 
4nm 2 to the minimum of 1 element per 0.5/zm 2 . Convergence, mesh independence, and 
accuracy was implemented using standard numerical procedures. Validation of the numerical 
results were established using homogeneous states. 



III. RESULTS AND DISCUSSION 



A. Linear regime: shape instability 



A sharp interface model for the isotropic/nematic mesophase front growth is presented 



in the Appendix. This derivation differs from standard sharp interface models 
that it: 

• accounts for different thermal conductivities of both phases. 
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2ij in 



there is no assumption that the ordered phase is isothermal 

uses an interfacial nematodynamic model [l5| for the velocity of the interface as a 
function of temperature (and not the typically used Gibbs-Thompson relation 241]). 



The final expression derived from the simple sharp-interface model (refer to the appendix) 
for the dispersion diagrams, in the low amplitude linear regime is: 
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where k = 2tt/\ is the wave vector, A is the wavelength, a is the growth coefficient, f n is 
the bulk nematic free energy at the interface temperature, (3 is a surface viscosity (defined 
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in the appendix), L is the heat of phase transition, a = k/pCp is the thermal diffusivity, 
v is the interface velocity, 7 is the surface tension. The superscripts n and % refer to the 
nematic and isotropic phases, respectively. 
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Figure 6: Sample order parameter and temperature profiles for initial flat interfaces from the 
sharp-interface model (solid lines) and stationary simulation results (dotted lines) a) (top) tem- 
perature versus y-coordinate over full domain b) (middle) uniaxial nematic scalar order parameter 
versus y-coordinate over full domain c) uniaxial nematic scalar order parameter versus y-coordinate 
magnified in the interface region. 

The sharp-interface model used to derive Eq. ffl5|) uses an important approximation: 
the value of the order parameter is assumed uniform in the nematic phase. This approxi- 
mation, although convenient for the derivation of an analytical solution, is not physically 
realistic. The order parameter in the bulk nematic phase changes due to its dependence on 
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temperature. This non-uniformity in the order parameter, observed in the simulations, has 
an effect on the temperature profile which is not accounted for in the sharp interface model. 
This effects results from the presence of temperature gradients which is accounted for in the 
energy balance (see term 3 of Eq. (jSJ)). Sample order parameter and temperature profiles 
are shown in Figure El for both the sharp-interface model and the non-isothermal Landau-de 
Gennes model to elucidate this point. 

Simulations were performed using different initial wavelengths, under two different tem- 
perature gradients, where the initial perturbation amplitude used was much less than the 
wavelength. The material properties and individual simulation parameters are given in 
Tables H] and [Til The boundary temperatures were selected so that the initially imposed 
interfacial temperature was equal for both gradients. The perturbation amplitude, defined 
as the difference of the vertical positions of the interface at the symmetry axes (see Fig. [5]) 
was obtained from the simulations as a function of time. The position of the interface was 
determined from the contour at which the uniaxial nematic scalar order parameter equaled 

Si = ^p- (16) 

where <S& bulk value of the uniaxial nematic order parameter that corresponds to the initial 
interface temperature (7$). The lack of biaxiality at the interface on the symmetry axes was 
verified a posteriori. Finally, the amplitude versus time was fitted with an exponential and 
the growth coefficient (a), was obtained: 

A(t) = A exp (at) (17) 

The results are shown in Fig. UJ, where the growth coefficient obtained from simulations 
for different wavelengths in the small-amplitude regime are compared with those obtained 
from the sharp-interface model using two different values of the surface energy: a value 
from previous work [40| and the values estimated from both sets of simulations. The sharp 
interface mode predicts a sharp decrease in the growth coefficient in the vicinity of k = 
which is shown in the insets of Fig. [7J 

As previously explained (see Sec. [TJ), a value of k exists for which the growth coefficient 
is maximum. For large values of k (low wavelength), curvature effects are predominant and 
thus the growth coefficient is negative (the perturbation shrinks). At intermediate values, the 
growth coefficient is positive because of the increased destabilizing effect of the temperature 
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Figure 7: Dispersion diagram results directly from simulation (circles), using the analytical solution 
(dashed line, calculated using surface energy obtained from the definition in ref 4.5 x lO -4 ^), 
and using the analytical solution (full line, dispersion diagram calculated using surface energy = 
3.0 x 10~ 5 -^ best fit from simulation). Insets in the upper right corner of the plots show initial 
regime (low n) are shown for reference, a) (top) results corresponding to simulations la — Id, b) 
(bottom) results corresponding to simulations 2a — 2d. 




Figure 8: (Color online) initially homeotropic transition from the linear regime i) t = 156.0//S ii) 
t = 175.5/ms iii) t = 195.0//S iv) t = 214.5/is; the surface corresponds to the scalar uniaxial nematic 
order parameter (red/black corresponds to isotropic/nematic), and the horizontal length scale of 
the figure is lfj,m. 
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Table I: Material properties for 5CB used in simulation 




PaTflTTipfpr 

A CI' J. CliLL±\J Uv^J- 




Units 




0.084 


N s /?n 2 




307.2 


K 




1.4 x 10 5 


J/m 3 K 


b 


1.8 x 10 7 


J/m 3 


c 


3.6 x 10 6 


J/m 3 


h 


3.0 x 10~ 12 


J 1 TYl 


h 


3.1 x 10~ 12 


J lm 


u 


0.0 x 10~ 12 


J/tti 


h 


0.2009 


W/mK 


k± 


0.1364 


W/mK 




1800 


J/kgK 


p 


1000 


kg/m 3 



Table II: Simulation parameters and growth coefficient (a) results. 



# 


B(fim) 


Ti(K) 


T 2 (K) 


X (fj,m) 


v(m/s) 


a(fjLm) 


Ti{K) 


^(xlO" 1 ^" 1 ) 


la 


15 


307.4 


305.0 


1 


0.01955 


0.06 


306.27 


-9.23 


lb 


15 


307.4 


305.0 


2 


0.0195 


0.06 


306.27 


4.134 


lc 


15 


307.4 


305.0 


4 


0.0195 


0.06 


306.27 


5.784 


Id 


15 


307.4 


305.0 


7.5 


0.0196 


0.09 


306.27 


5.282 


2a 


15 


306.8 


305.4 


1 


0.0202 


0.06 


306.23 


-9.744 


2b 


15 


306.8 


305.4 


2 


0.02035 


0.06 


306.23 


3.467 


2c 


15 


306.8 


305.4 


4 


0.02045 


0.06 


306.23 


4.979 


2d 


15 


306.8 


305.4 


7.5 


0.02005 


0.09 


306.23 


4.615 



gradient at the interface. As the value of k approaches zero, when the wavelength becomes 
much larger than the characteristic length of the temperature profile, the interface behaves 
as a flat interface and the problem becomes essentially one-dimensional. When the external 
length B (see Fig. [5]) is larger than the characteristic thermal length (as in the present case), 
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the temperature at the interface in this one-dimensional problem depends only slightly on 
the position of the interface. Thus, recalling Eq. [TJ the growth coefficient will be small. 
As B approaches oo (as usually assumed in theoretical studies), the temperature at the 
interface becomes independent of position and the growth coefficient is zero. 

The parameters involved in the calculation of the sharp-interface dispersion diagrams have 
been calculated from simulation results, as in ref |40| (see Fig. [7]). The calculated dispersion 



diagram with these parameters reproduces well the maximum value of the growth coefficient, 
as compared with simulations, but the curve resulting from the surface tension definition 

n 

from ref 40] is more narrow than that determined from direct numerical simulation. 



B. Non-linear regime: structural dynamics 

The shape instability due to the presence of temperature gradients was not the only 
interfacial instability that was observed. A process of defect formation, shedding, and the 
growth of a disoriented planar domain was also observed following the initial linear regime. 
It is important to note that results presented in the previous section were in the transient 
linear instability regime. The nonlinear regime instabilities presented in this section were 
observed to follow the initially linear instabilities 

These nonlinear regime instabilities are shown in Figs. IBlfTUl which correspond to simula- 
tion 2c (see Table HT|) . Referring to Fig. [SJ a cusp is formed at the interface as the instability 
transitions from the linear regime, discussed in the Sec. IIII Al to the nonlinear regime. As 
the cusp sharpens, it forms a +| disclination which is then shed into the bulk, similar to 
the defect shedding mechanism found by Wincure and Rey in growing initially homogeneous 
nematic spherulites under isothermal conditions [3] . As the defect sheds, a planar anchored 
interfacial regime nucleates and grows with the moving front. Fig. [TOa shows the full view of 
the computational domain with the director profile. Fig. [TUb shows the temperature profile 



and biaxia 
as follows 



nature of the +5 disclination and planar isotropic/nematic interface, computed 



44|: 



P 2 = l-6 l ^^^7 l (18) 



[( Q Q) Q ] 

(Q : Qf 

ranging from ranging from fully uniaxial f3 2 = to fully biaxial (3 2 = 1. 

Biaxiality has a key role in this process, where the transition to an inherently biaxial 
planar interface from a uniaxial homeotropic interface is achieved through the shedding of 
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a biaxial disclination defect. Exhaustive past work has addressed these interfacial texturing 
processes for isothermal isotropic/nematic interfaces, see refs. jlol . 40 ] for a full treatment. 
This defect formation, shedding, and planar interfacial growth phenomena are driven by 
planar anchoring having a lower energy than that of a homeotropic anchoring. This is 
shown in Fig. [TTJ which shows the gradient free energy density f g (terms 4-5 Eq. ([3])): 

(19) 



fg = ^i(VQiVQ) + \h (V ■ Q) ■ (V ■ Q) 



across the interface at the horizontal position at which the defect forms, as a function of 
the position (across the interface) and time. Before the defect formation, the interfacial 
anchoring is homogeneous in the high-energy homeotropic state. As the disclination forms, 
the interfacial anchoring becomes heterogeneous where an energetically favorable planar 
regime forms trailing the shed defect. In the vicinity of the disclination defect the free 
energy is greatly increased; it is important to note that the total free energy of the system 
is minimized in that the energetically favorable interfacial region of planar anchoring grows, 
decreasing the total free energy. 



( j ) 1 i I i 1 1 i i i i i i i i i i i i 
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Figure 9: (Color online) initially homeotropic defect formation i) t = 214.5/is ii) t = 219.4/is iii) 
t = 229. l/us; the surface corresponds to the scalar uniaxial nematic order parameter (red/black 
corresponds to isotropic/nematic), arrows correspond to the uniaxial nematic director (should be 
considered headless), and horizontal length scale is 0.65/um 



In order to exhaustively confirm that this instability results from the difference in an- 
choring energies, a simulation was performed with Z 2 = 0. This is based upon the notion 
that for an interface, neglecting curvature and biaxiality, the surface energy is predicted to 



be 



45|: 



7 



b 3 ^3h + l 2 /2 + 3l 2 (n-k) 2 /2 



(20) 



486c 5 / 3 

where k is the unit vector normal to the interface. Equation (|20|) shows that the sign of 
l 2 determines which orientation has the lower surface energy; if it is positive it will planar 
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(n _L k), and if it is negative, homeotropic (n = k). When I2 = 0, neither of the interfacial 
anchorings are preferred. Thus if the surface free energy is the driving force for defect 
formation; no defect would be expected to form when l 2 = which was confirmed via this 
simulation where, for long times, no defect formation was observed. 




Figure 10: (Color online) t = 229.1/xs (corresponds to Fig. [S](iii) a) (left) view of full computational 
domain where the surface corresponds to the scalar uniaxial nematic order parameter (red/black 
corresponds to isotropic/nematic), arrows correspond to the uniaxial nematic director (should be 
considered headless), and the horizontal length scale is 2fim b) (right) magnified view of the planar 
interface with the surface corresponding to j3 2 (see Eq. (|18j) ) and the contours corresponding to 
temperature (the minimum/maximum/increment is 306. 285-ftT/306. 2951^/0.002^) the horizontal 
length scale is 0.6/im. 



IV. CONCLUSIONS 

In this work, two-dimensional multi-scale simulation of the isotropic/nematic mesophase 
was used to study both linear and nonlinear morphological instabilities under directional 
growth. The major contributions of this work are: 

• an extended sharp-interface model was derived (refer to Sec. IIII Al and the appendix) 
and applied to linear morphological instabilities of the isotropic/nematic mesophase 
transition. 

• results from comparing the sharp-interface model and the two-dimensional simulation 
of the tensorial Landau-de Gennes model (refer to Sec. Ill Bj) identify the limitations of 
sharp-interface models to predict the maximum value of the growth coefficient (see Fig. 

21 




Figure 11: (Color online) Surface plot of the gradient free energy density (refer to Eq. ([3])) 
as a function of time and y-coordinate (jum); features of the energy surface include the initial 
homeotropically oriented interface, defect shedding (the high energy peak), and the relatively 
lower energy planar oriented interface. 

[7j). It is shown that this type of model fails to adequately capture the time-dependent 
morphological instability growth compared to full two-dimensional simulation of a 



non-isothermal Landau-de Gennes model 35]. 



nonlinear instabilities were studied for an initially homeotropic isotropic/nematic in- 
terface, where the phenomenological model predicts the planar isotropic/nematic in- 
terface being more stable. Disclination shedding at the interface, similar to experi- 

n n 

mentally [22[ and theoretically observed phenomena [16], and the formation/growth 



of bulk texturing were found to result (see Figs. ISl fTTj) . 

While these results confirm the validity of sharp-interface models for the well-studied prob- 
lem of linear morphological instabilities, a substantial conclusion can be drawn from simu- 
lation results in the nonlinear regime. That conclusion being that, using the full tensorial 
Landau-de Gennes model for the isotropic/nematic phase transition, extended to the non- 
isothermal case [35J], nonlinear morphological instabilities can be accessed that are both 
inaccessible using more simple theoretical approaches and difficult to characterize experi- 
mentally. 
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Furthermore, three-dimensional simulation using the tensorial Landau-de Gennes model 



351 ] . based upon the current results, will enable the study of geometries that are more di- 



rectly representative of experimentally conditions. For example, some of the physics that 
can be accessed include surface anchoring, meniscus formation [17J, three-dimensional tex- 
tures/defects et al. 
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Appendix A: SHARP INTERFACE MODEL FOR LINEAR STABILITY OF A 
PLANAR INTERFACE GROWING IN AN EXTERNAL TEMPERATURE GRA- 
DIENT 

The objective of this appendix is to derive an analytical expression for the dispersion 
diagram, considering a sharp-interface model. The energy balance for the bulk nematic and 
isotropic phases, in a frame moving with velocity v are: 

dT n Qrpn /Q2 T n Q2 T n\ 

V ^T = ® " -5T5- + -5I5- ( Al ) 



dt dy \ dy 2 dx 



dT i dT i , d 2 T i d 2 T i~ 

*-*r = <* -SIT + -5-5- ( A2 ) 



dt dy \ dx 2 dy 2 

where a is the thermal diffusivity. The superscripts % and n refer to the isotropic and 
nematic phase, respectively. We are assuming that the order parameter is uniform in the 
whole nematic phase. 

The boundary conditions are: 

T\ y= „ B = Ti T\ y=B = T 2 (A3) 

and, additionally, at the interface y = h(x) the temperature must be continuous and the 
heat released by the phase transition must equal the net heat flux. 

For a steady-state flat interface (base solution), we take the velocity of the moving ref- 
erence system v as the velocity of the interface, and place the interface at the origin, so 

23 



h(x) = 0. The temperature is only a function of y, and the solution to Eqs. (lAlj) and (1A2|) 
is: 



T l = C* jexp + Tj (A4) 

where C n and C l are integration constants, and T/ is the temperature at the interface. The 
boundary conditions Eq. IA3I and previously mentioned interface conditions are then used 
to find the particular solution from the general solution Eq. IA4[ 

Tj = {T x (/(-«) - 1) - T 2 (/(«)- 1) 

+ ^ (/(«) - 1) (/HO - !)}/{/(-«) " /(«)} 

= exp f —5 J ; f n (v) = exp ( —5 



Oii I \OL 



'II 



C n = H II C i = 12 II , A5) 

exp(^)-l' eX p(-X)_i 1 } 
Now a perturbation to the base solution is considered: 

rpi rpi _|_ rp*i 

T n = T n + T *n ^ 

where the T* 1 and T* n are the perturbations, and the subscript b are the base solutions. 
Replacing Eq. (jXfiJ in Eqs. ((AT]) and flA2l) : 

9T *n 5T *n fd 2 T* n d 2 T* n \ /A , 

«-5T- = « "HJ" + -53- ( A7 ) 



<9t ch/ \ dy 2 dx 



dt dy a \ dy 2 dx 2 ) ^ ^ 

As long the perturbations are small and do not reach the borders of the simulation domain, 
boundary conditions at infinity can be used: 

T* n \ =0 T H \ =0 (A9) 

>y— — oo \y=oo v ' 

Writing the boundary conditions at the interface, considering that now h = h(x) and lin- 
earizing, the following boundary conditions are obtained, where the temperatures and their 

24 



derivatives are evaluated at y = 0: 



r)T n 
dy 

+ 



dT[ 
dy 



dy d 2 y 



h + T* 1 
dy ~^ d 2 y 



(A10) 



dh 



(All) 



The velocity of the nematic-isotropic interface, can be calculated from a nematodynamic 
interfacial model 151 ]: 

dQ s 



(3w=(L + V s -T s )-k + ^Q s : 

at 

= /a„ : dA 



(A12) 
(A13) 



where (3 is the interfacial viscosity, L is the temperature dependent net stress load (which 
reduces to the free energy difference between the nematic and the isotropic phase), V s -T s is 
the capillary force (neglecting anisotropy in surface tension, it is the product of the surface 
tension and the bidimensional curvature of the interface), and the last term is the change 
in the value of the order parameter at the interface (assumed to be 0). If the perturbation 
is small, Eq. IA12I reduces to: 



(3w = -r + 



-fd 2 h 

2~cb? 



(A14) 



where f n is the bulk nematic free energy, and 7 is the surface tension. The surface tension, 
neglecting curvature, is found to be: 



7 



df n 
dVQ 



: (VQ) T d\ 



(A15) 



where Ig is the surface 2x2 identity matrix. After linearization, considering that for small 
perturbations the normal to the interface is approximately the y direction, the change in the 
velocity with respect to the base solution (the velocity of the interface in the moving frame) 
is: 



dh 
It 



dT 



<9T™ 



dy 



h + T* r 



y=o 



\y=0 



+ 



7 d 2 h 
2pd^ 



(A16) 



where the derivative with respect to temperature is a total derivative and thus the equilib- 
rium order parameter as a function of the temperature must be used in the free energy. 
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The solution, for a sinusoidal perturbation at the interface, is: 



h 



replacing in Eq. flA7j) results in: 



a 



a 



Tq exp(aj?/) exp(crt) sm(K,x) 
Tq exp(a n y) exp(cxi) sin(/tx) 
hoe^ sm(nx) 

— VCLi = Oii [o^ — K 2 ) 

- vai = a n (a 2 n - k 2 ) 



(A17) 



(A18) 



The system formed by Eqs. (1A10[ IA161 IA17[ IA18I) is a homogeneous system. In order to 
have a non trivial solution, the following dispersion relation must be satisfied: 







no-l\ 



d(f n P 



dT 



a C n + — 

a 2p 



d(f n (3 
dT 



--U n ( — ) C n + ai (-) CntSn + S,)- 1 

o \a n J \atij J 



(A19) 



+ 2 o- 



where: 



\jv 2 + 4a j (k 2 + a), j = i,n 



(A20) 
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